import numpy as np
import matplotlib.pyplot as plt
Question 1. — Créer une matrice de taille 8 8 contenant des 0 et des 1 à la manière d’un échiquier.
Z = np.zeros((8,8))
Z[1::2,::2] = 1
Z[::2,1::2] = 1
Z
Question 2. — Définir la matrice M = 1 5 9 13 17 2 6 10 14 18 3 7 11 15 19 4 8 12 16 20 1
et en extraire la matrice 0
6 18 10 7 19 11 5 17 9 1
M = np.arange(1,21)
M = M.reshape((5,4))
M = M.transpose()
print(M)
print(M[1:,:][:,[1,4,2]])
On se place dans un cadre de régression en dimension 1, c’est-à-dire l’ensemble des entrées est X = R et l’ensemble des sorties Y = R. Les prédicteurs sont donc des fonctions de la forme f : R => R.
Dans ce TP, nous n’allons pas travailler avec des données réelles, mais nous allons les générer nous-mêmes. On considère la fonction polynomiale
$g(x) = \frac{3}{2} x^3 - x^2 - \frac{3}{4} x + 1 $
qu'on va utiliser pour générer un échantillon d'apprentissage $(X_i, Y_i)_{i\in[1,m]}$ i.i.d selon :
$X_i \sim Unif([-1,1])$
$Y_i = g(X_i) + \frac{1}{20} \zeta_i$
où $\zeta_i \sim \mathcal N(0,1)$
Question 3. — Définir un array x contenant m = 15 réels tirés uniformément sur [-1; 1]. On pourra utiliser la fonction np.random.random_sample.
m = 15
x = 2*np.random.random_sample(m)-1
x
Question 4. — Définir un array y contenant les valeurs Yi i appartenant à n définies selon (2.2). On pourra faire appel à la fonction np.random.randn.
def g(x):
return 1.5*x**3 - x**2 - .75*x + 1
ksi = np.random.randn(m) #suit une n(0,1)
y = g(x) + .05*ksi
y
Question 5. — Visualiser dans le plan les données obtenues (Xi; Yi) i appartenant de 1 à m sous la forme de points bleus.
plt.plot(x,y,'bo')
plt.show()
Question 6. — Avec le même procédé, générer un échantillon de test (Xi',Yi') i appartenant à 1:m' de taille m' = 30 qu’on stockera dans des array x_test et y_test.
m_test = 30
x_test = 2*np.random.random_sample(m_test)-1
y_test = g(x_test) + .05*np.random.randn(m_test)
y_test
Nous allons faire appel à la bibliothèque scikit-learn pour construire des prédicteurs à partir des données d’apprentissage. Dans scikit-learn, tout prédicteur possède une fonction fit qui permet de spécifier les données d’apprentissage, ainsi qu’une fonction predict qui permet, une fois le prédicteur entraîné, de calculer la prédiction de le prédicteur sur de nouvelles entrées. On cherche ici à construire le prédicteur linéaire (c’est-à-dire affine) qui minimise l’erreur des moindres carrés. Autrement dit, la classe de prédicteurs considérée est
$\mathcal F = \{f_{a,b} : x \longmapsto ax+b\}_{a,b \in R}$
et on cherche à calculer le minimiseur du risque empirique pour la fonction de perte :
$l(y,y') = (y-y')^2$
$$\hat f = \underset{f \in \mathcal F}{argmin} \{ \sum_{i=1}^m (Y_i - f(X_i))^2 \} $$Nous allons pour cela utiliser le module LinearRegression que l’on charge comme suit.
from sklearn.linear_model import LinearRegression
f = LinearRegression()
On a initialisé ci-dessus un prédicteur f de régression linéaire aux moindres carrés. Il faut à présent fournir les données d’apprentissage via la fonction f.fit. Celle-ci exige que les données d’entrée ($Xi)_{i \in m}$ soient présentés sous la forme d’un array de dimension 2 contenant m lignes (une par donnée) et 1 colonne (car ici l’espace d’entrée X = R est de dimension 1). Or, x comporte une seule dimension. Nous allons donc définir X de sorte qu’il contienne les valeurs de x sous la bonne forme, et de même pour x_test
En somme, Le modèle s'attend à ce qu'on lui donne :
ici, comme X est dans R et non dans R^2, le modèle ne va pas comprendre donc on ajoute une dimension en plus
x
X = x[:,np.newaxis]
X
#le np.newaxis permet d'ajouter une dimension en plus
#comme le modèle s'attend à du 2 dimensions et que X est juste un vecteur (donc une dimension (m,) ),
#nous utilisons np.newaxis pour ajouter une dimension (m,1)
#Exemple : X = [1, 2, 3, 4]
#1) X[:,np.newaxis] "transforme" ton X en dimension (m,1)
#[ [1, 2, 3, 4] ]
#2) X[np.newaxis,:] "transforme" ton X en dimension (1,m)
#[ [1],
#[2],
#[3],
#[4] ]
#X.reshape((len(X), 1)
#attention, dans notre cas : "dimension" = (num ligne, num colonne) de notre tableau
print(x) #15 x 1
print(x.shape)
print(X) #2 x 15
print(X.shape)
X_test = x_test[:,np.newaxis]
X_test
On peut à présent spécifier les données d'apprentissage.
f.fit(X,y)
Si on écrit $\hat f(x) = \hat ax + \hat b$, les coefficients $\hat a$ et $\hat b$ du prédicteur f sont donnés par f.coef et f.intercept respectivement. Les prédictions du prédicteur f sur les données d’apprentissage d’une part, et données de test d’autre part sont donnés par:
f.predict(X)
f.predict(X_test)
Question 7. — Produire une figure qui superpose aux points (Xi; Yi) i appartenant à 1:m la droite correspondant au graphe du prédicteur $\hat f$.
xplot = np.linspace(-1,1,500).reshape(-1,1)
plt.plot(x,y,'bo') #les points
plt.plot(xplot,f.predict(xplot)) #la droite
plt.show()
Remarque : la constante est un paramètre que ton modèle détermine "tout seul" en minimisant le risque la méthode "fit" permet d'entraîner ton modèle d'apprentissage grâce aux données d'entraînement données en argument
#Même chose sans constante, on n'obtient pas la même chose ?
g = LinearRegression(fit_intercept=False)
g.fit(X,y)
xplot = np.linspace(-1,1,500).reshape(-1,1)
plt.plot(x,y,'bo')
plt.plot(xplot,g.predict(xplot))
plt.show()
Question 8. — Calculer l’erreur moyenne d’apprentissage ainsi que l’erreur moyenne de test. Commenter.
print(sum((y-f.predict(X))**2)/m)
print(sum((y_test-f.predict(X_test))**2)/m_test)
L'erreur de test est plus importante que celle d'apprentissage ?

On souhaite dans cette question calculer le prédicteur $\hat f$ pour $n = 2$. Nous allons implémenter la fonction $\psi$ introduite ci-dessus à l’aide d’outils fournis par scikit-learn.
from sklearn.preprocessing import PolynomialFeatures
psi = PolynomialFeatures(2,include_bias=False).fit_transform
# include_bias = pas de constante
#"fit_transform" est une fonction (qui est vide pour l'instant) afin de fit les
# données qu'on lui donne plus tard en paramètres
#psi = .... .fit_transform # on a définit une fonction
# psi(X) # on applique cette fonction sur X
#comme dans LinearRegression vous avez déjà inclus le biais (intercept)
# donc pas besoin de l'ajouter aussi dans PolynomialFeatures
#dans l'énoncé, il est dit que psi(x) vaut [x, x^2, ...],
# nous n'avons pas [1, x, x^2, ..] car on a mis
#include_bias=False
PolynomialFeatures : Generate polynomial and interaction features. Generate a new feature matrix consisting of all polynomial combinations of the features with degree less than or equal to the specified degree. For example, if an input sample is two dimensional and of the form [a, b], the degree-2 polynomial features are [1, a, b, a^2, ab, b^2].
option : include_bias : boolean. If True (default), then include a bias column, the feature in which all polynomial powers are zero (i.e. a column of ones - acts as an intercept term in a linear model).
f = LinearRegression().fit(psi(X),y)
plt.plot(x,y,'bo')
plt.plot(xplot,f.predict(psi(xplot)))
plt.show()
print(sum((y-f.predict(psi(X)))**2)/m)
print(sum((y_test-f.predict(psi(X_test)))**2)/m_test)
On observe une nette amélioration par rapport à la regression linéaire
Question 10. — On souhaite à présent généraliser le calcul à tout n > 1 (qui correspond au degré maximal des polynômes).
a) En s’inspirant de la question précédente, écrire une fonction qui prend en argument l’entier n, et qui renvoie deux éléments: le prédicteur $\hat f$ et la fonction $\phi$.
# Question a
def reg_poly(n):
psi = PolynomialFeatures(n,include_bias=False).fit_transform
return LinearRegression().fit(psi(X),y), psi
b) Calculer le prédicteur pour n appartient à {3, 4, 13, 14} et le visualiser (on pourra si besoin spécifier des limites pour l’abscisse).
#Question b
poly = {}
psi = {}
plt.plot(x,y,'bo')
for n in [3,4,13,14]:
poly[n], psi[n] = reg_poly(n)
plt.plot(xplot,poly[n].predict(psi[n](xplot)))
plt.axis([-1, 1, -1.5, 1.5])
plt.legend(['données','n = 3', 'n = 4', 'n = 13', 'n = 14'])
plt.show()
c) Pour n = 3, comparer les coefficients du prédicteur avec ceux de la fonction g qui a généré les données.
# Question c
print(poly[3].intercept_)
print(poly[3].coef_)
Pour n=3, ces coefficients sont proches à ceux de la fonction g
Rappel : $g(x) = 1,5 x^3 - x^2 - 0,75 x + 1 $
d) Que se passe-t-il pour n = 14 ?
#Question d
plt.plot(x,y,'bo')
for n in [14]:
poly[n], psi[n] = reg_poly(n)
plt.plot(xplot,poly[n].predict(psi[n](xplot)))
plt.show()
Pour n=14, il existe un polynôme qui passe par chacun des 15 points (polynôme d'interpolation de Lagrange). Ce polynôme minimise évidemment le risque empirique, c'est donc celui donné par la minimisation du risque empirique
e) Tracer une figure avec en abscisse n = 1;... 14, et en ordonnée les graphes des erreurs moyennes d’apprentissage et de test. Commenter.
#Question d
mean_train_error = [] #erreur d'apprentissage
mean_test_error = [] #erreur de test
n_max = 14
for n in range(1,n_max+1): #n va de 1 à 14 (abscisse)
f, psi = reg_poly(n)
mean_train_error.append(sum((f.predict(psi(X))-y)**2)/m)
mean_test_error.append(sum((f.predict(psi(X_test))-y_test)**2)/m_test)
plt.plot(range(1,n_max+1),mean_train_error)
plt.plot(range(1,n_max+1),mean_test_error)
plt.axis([0, n_max, 0, .2])
plt.show()
Plus n est élevé :
La grande complexité des polynôme permet de mieux coller aux points d'apprentissage, mais ne permet pas une bonne généralisation à de nouveaux points.
On considère à présent la régression polynomiale aux moindres carrés avec régularisation LASSO:
$$\hat f = \underset{f \in \mathcal F_n^\text{poly}}{argmin} \{ \frac{1}{2m} \sum_{i=1}^m (Y_i - f(X_i))^2 + \alpha \sum_{k=1}^n | a_k | \} $$où $\alpha$ > 0 est un paramètre à choisir.
Avec sklearn, on initialise une régression linéaire avec régularisation LASSO de la façon suivante où l'arguement correspond au choix du paramètre $\alpha$ :
from sklearn.linear_model import Lasso
f = Lasso(0.1) #alpha = 0.1
Question 11. En s’inspirant de la section précédente, calculer pour n = 14, la régression polynomiale avec régularisation LASSO pour différentes valeurs du paramètre (par exemple $\alpha \in \{ 10^{-6}, 10^{-5}, \dots , 10^{-1}, 1 \}$), et les visualiser. Que remarque-t-on sur les polynômes obtenus ? sur leurs coefficients ? leurs degrés ?
from sklearn.linear_model import Lasso
def reg_poly_lasso(n, alpha=1):
psi = PolynomialFeatures(n,include_bias=False).fit_transform
return Lasso(alpha, max_iter = 1000000).fit(psi(X),y), psi #max iter new
for alpha in [10**a for a in range(-6,1)]:
f, psi = reg_poly_lasso(14,alpha)
plt.plot(x,y,'bo')
plt.plot(xplot,f.predict(psi(xplot)))
plt.axis([-1, 1, -1.5, 1.5])
print('Coefficients:',f.coef_)
plt.show()
Pour des valeurs de alpha assez grands, l'algorithme produit des polynômes ayant beaucoup de coefficients nuls et de faible degré. Cela limite la complexité des polynômes et donc le sur-apprentissage.
Question 12. — Tracer les erreurs d’apprentissage et de test pour les différentes valeurs de $\alpha$ (on pourra utiliser une échelle logarithmique en abscisse afin d’améliorer la lisibilité). Commenter.
mean_train_error = []
mean_test_error = []
nonzero_coefficients = []
degree = []
alpha_range = [10**a for a in range(-6,1)]
for alpha in alpha_range:
f, psi = reg_poly_lasso(14,alpha)
mean_train_error.append(sum((f.predict(psi(X))-y)**2)/m)
mean_test_error.append(sum((f.predict(psi(X_test))-y_test)**2)/m_test)
print("erreurs d'entraînement : ", mean_train_error,"\n\n")
print("erreurs de test : ", mean_test_error)
plt.plot(np.log10(alpha_range),mean_train_error)
plt.plot(np.log10(alpha_range),mean_test_error)
plt.axis([-6,0,0,0.5])
plt.legend(['erreur de train','erreur de test'])
plt.show()
Le paramètre alpha joue un rôle similaire à n dans la question 10 (mais en sens inverse).
Lorsque alpha décroît, les polynômes obtenus sont de plus en plus complexes, et l'erreur d'apprentissage diminue. Mais quand alpha est assez petit, l'erreur de test augmente à nouveau: il y a alors sur-apprentissage.
Nous allons travailler avec le jeu de données de Ronald Fischer, qui date de 1936, et qui regroupe des données sur des iris.
import numpy as np
import seaborn as sns
import pandas as pd
import matplotlib.pyplot as plt
iris = sns.load_dataset("iris")
iris contient alors les données sous la forme d’une DataFrame: il s’agit d’un type de données fourni par la librairie Pandas, qui permet de stocker des données sous la forme d’un tableau dont les colonnes sont nommées, et qui fournit un grand nombre de fonctions pour explorer et manipuler ces données. On peut afficher les 5 premières lignes avec la commande suivante.
iris.head()
On voit que les 4 premières colonnes correspondent respectivement à la longueur du sépale, la largeur du sépale, la longueur de la pétale et la largeur de la pétale (toutes ces mesures sont en centimètres). Ces 4 variables vont constituer l’entrée. La cinquième colonne correspond à l’espèce de l’iris (qui peut être: Setosa, Versicolor ou Virginica) et constitue la sortie. Il s’agit donc d’un problème de classification: on souhaite construire un prédicteur prédisant l’espèce en fonction des autres caractéristiques. Par ailleurs, la DataFrame comporte 150 lignes: il s’agit du nombre de données.
Nous souhaitons dans un premier temps simplifier le jeu de données en réduisant la dimension de l’espace des entrées de 4 à 2. Il s’agit donc de déterminer quelles sont les deux variables explicatives qui semblent les plus prometteuses pour la prédiction de l’espèce. Pour ce faire, nous allons utiliser la fonction sns.pairplot de la librairie Seaborn.
sns.set()
sns.pairplot(iris, hue="species");
plt.show()
Question 13. — Choisir les deux variables explicatives qui semblent les plus prometteuses pour la prédiction de l’espèce, c’est-à-dire celles qui dans la figure semblent le mieux séparer les différentes espèces. Créer un array X contenant les colonnes correspondant aux deux variables explicatives sélectionnées, ainsi qu’un array y contenant la colonne correspondant à l’espèce (on pourra faire appel à iris.values qui renvoie l’ensemble des données de la DataFrame iris sous la forme d’un array NumPy).
=> On peut choisir un peu ce qu'on veut. On aimerait que les fleurs de la même espèce soient voisins. Quand on fait le graph de petal_with par rapport à petal_length ou a une super séparation. On choisit celles qui semblent le mieux différencier les trois classes sur le pairplots donc très corrélées à la variable "classe" en qq sorte.
Argument du chargé de TD Hu (faux ?) Les variables sepal_width et sepal_length semblent très corrélés (idem pour les deux variables petal). Donc l'idée est de sélectionner une variable sepal et une variable petal
Les deux variables explicatives qui semblent le mieux séparer les espèces sont petal_width et petal_length
#sepal_length sepal_width petal_length petal_width species
X = iris.values[:,2:4] #petal_length, petal_width Attention : 2:4 = 2 et 3 !
y = iris.values[:,4] #species
Nous souhaitons à présent séparer le jeu de données en deux échantillons: un d’apprentissage de taille 90, et un de test de taille 60. Une première idée serait d’utiliser les 90 premières données pour constituer l’échantillon d’apprentissage et les 60 dernières pour l’échantillon d’entraînement.
X_train = X[:90]
y_train = y[:90]
X_test = X[90:]
y_test = y[90:]
On entraîne le classifieur sur les données d'apprentissage et on teste sur jeu de test.
Warning : on ne teste pas le prédicteur sur le jeu d'entraînement car on minimise l'erreur sur le jeu d'entraînement donc ce sera forcément biaisé !
Question 14. — Observer les échantillons ainsi obtenus et expliquer en quoi ils sont problématiques. On peut remédier au problème en faisant un shuffle
Les données X et y n'étaient pas mélangées: les espèces sont regroupées. En particulier, l'échantillon de test ne contient quasiment que des espèces virginica, alors que l'échantillon d'apprentissage ne contient aucune espèce virginica. Cela fausserait donc et l'apprentissage et le test. Il faut donc au préalable mélanger le jeu de données.
Warning : on suppose que distribution de jeu d'entraînement et jeu de test sont les mêmes. C'est pourquoi il est important de mélanger les données avant de les séparer. Si on n'avait pas mélangé, le prédicteur n'aurait vu que des setosa et n'aurait pas plus prédire autre chose.
from sklearn.utils import shuffle
X, y = shuffle(X,y)
X_train = X[:90]
y_train = y[:90]
X_test = X[90:]
y_test = y[90:]
On peut à présent entraîner un prédicteur kNN, pour k = 3 par exemple.
La distance utilisée par défaut est la distance euclidienne. Les prédicteurs de classification de scikit-learn fournissent une fonction .score (donc ici knn.score) qui donne la proportion de bonnes prédictions sur un échantillon fourni en argument.
from sklearn.neighbors import KNeighborsClassifier
knn = KNeighborsClassifier(n_neighbors=3)
knn.fit(X_train, y_train)
print("fraction de prédictions bonnes (train) : ",knn.score(X_train,y_train))
print("fraction de prédictions bonnes (test) : ",knn.score(X_test,y_test))
Question 15. — Quelle est la relation entre le score défini ci-dessus et le risque empirique?
Le score empirique mesure la proportion de bonnes prédictions, alors que le risque empirique mesure la proportion de mauvaises prédictions (lorsqu la fonction de perte considérée est la perte 0-1). Donc: score = 1 - risque empirique.
Le risque : proba de se tromper dans la prédiction et donc empiriquement c'est la moyenne de l'indicatrice de on s'est trompé.
Le score : le nombre de fois où on a eu raison.
Question 16. — Calculer le prédicteur kNN pour $k \in \{1, \dots, 20\}$. Tracer les courbes des scores calculés sur les échantillons d’apprentissage et de test respectivement. En déduire une valeur de k qui semble avoir produit le meilleur prédicteur.
knn = {}
score_train = []
score_test = []
k_range = range(1,21)
for k in k_range:
knn[k] = KNeighborsClassifier(n_neighbors=k)
knn[k].fit(X_train,y_train)
score_train.append(knn[k].score(X_train,y_train))
score_test.append(knn[k].score(X_test,y_test))
plt.plot(list(k_range),score_train)
plt.plot(list(k_range),score_test)
plt.legend(['erreur de train','erreur de test'])
plt.show()
Version Hu : La valeur k=6 semble ici donner le meilleur score de test.
Version Solène : On remarque que le score (bleu) reste relativement bon sur le jeu d'entrainement. Sur le jeu de test (orange), il atteind un plateau pour k entre 3 et 15. On choisit la valeur k=6
(Dernière version Kim) : Ici on choisirait plutôt k minimal...
Question 17. — En reprenant le code écrit précédemment, définir une fonction best_knn_score qui prend en argument deux array X et y, construit les échantillons d’apprentissage et de test, calcule les prédicteurs kNN pour $k \in \{1, \dots, 20\}$ et qui finalement renvoie le meilleur score calculé sur l’échantillon de test. Exécuter plusieurs fois cette fonction. Que remarque-t-on ? À quoi cela est-ce dû ?
def best_knn_score(X,y):
k_range=range(1,21)
score_test = []
X, y = shuffle(X,y)
X_train = X[:90]
y_train = y[:90]
X_test = X[90:]
y_test = y[90:]
for k in k_range:
knn = KNeighborsClassifier(n_neighbors=k)
knn.fit(X_train,y_train)
score_test.append(knn.score(X_test,y_test))
return max(score_test)
print(best_knn_score(X,y))
print(best_knn_score(X,y))
print(best_knn_score(X,y))
Le résultat est différent à chaque fois. Cela est dû à au mélange aléatoire des données par shuffle(X,y).
Le meilleure score sur l'échantillon de test varie, car le découpage des données entre test et apprentissage est aléatoire.
Question 18. — Écrire une fonction best_knn_score_avg qui prend en argument des array X et y, qui exécute 100 fois best_knn_score(X,y) et qui renvoie la moyenne des 100 scores obtenus.
def best_knn_score_avg(X,y):
scores = []
for i in range(100):
scores.append(best_knn_score(X,y))
return np.mean(scores)
print(best_knn_score_avg(X,y))
L’algorithme kNN est sensible à l’échelle utilisée pour chaque variable explicative. Pour s’en rendre compte, nous allons modifier l’échelle d’une variable explicative et observer la conséquence sur la qualité du prédicteur obtenu. On commence par créer une copie X de X. X = X.copy() La commande X = X n’aurait pas créé un copie, mais simplement un second nom pour le même objet: les modifications sur X auraient alors également affecté X.
Question 19. — Modifier l’array X en convertissant les données d’une des deux variables explicatives des centimètres aux mètres. L’utilisation de X au lieu de X affecte-t-elle la qualité des prédicteurs construits? Recommencer en convertissant cette fois-ci l’autre variable explicative en mètres (et en gardant la première en centimètres). Lorsqu’on souhaite s’assurer que l’influence d’une variable explicative ne soit ni trop faible, ni trop élevée en raison de son échelle, il est recommandée de normaliser les données. Autrement dit, pour chaque variable explicative, il s’agit d’appliquer une transformation affine telle que les valeurs qui en résultent aient une moyenne nulle et un écart-type égal à 1. scikit-learn propose une fonction pour cela. from sklearn import preprocessing X_scaled = preprocessing.scale(X)
X_ = X.copy()
X_[:,0] = X_[:,0]/100 #on change juste la première variable
print(best_knn_score_avg(X_,y))
Le score reste similaire
#X_ = X.copy()
X_[:,1] = X_[:,1]/100 #on change la 2e variable
print(best_knn_score_avg(X_,y))
Le score est cette fois-ci détérioré. (pas chez moi...)
On remarque que le fait de modifier l'échelle d'une variable explicative modifie les perfomances du prédicteur. C'est logique, puisque les plus proches voisins sont basés sur la distance euclidienne.
Question 20. — La normalisation des données donne-t-elle ici un meilleur résultat?
from sklearn import preprocessing
X_scaled = preprocessing.scale(X)
print(best_knn_score_avg(X_scaled,y))
Normaliser : variance des variables explicatives est la même.
En l'occurence, cela n'améliore pas le score, il est même légèrement détérioré. Oui, le score du predicteur est meilleur en moyenne si les données sont normalisées. => En fait, ça dépend des fois !
Question 21. — Si on considère l’ensemble des 4 variables explicatives initialement disponibles, le prédicteur obtenu est-il meilleur?
X = iris.values[:,0:4]
y = iris.values[:,4]
print(best_knn_score_avg(X,y))
Le score est meilleur. Le score n'est pas meilleur. => Ca dépend des fois en fait, donc non pas toujours ?